using System;

namespace NumericalRecipes
{
	/// <summary>
	/// Taken from Numerical Recipes in C++ and customized.
	/// </summary>
	public class Nr
	{
		public Nr(){}
		//	Static memory for speed
		static CholFloat statCholFloat=new CholFloat();
		static EigsDouble statEigsDouble=new EigsDouble();
		static EigsFloat statEigsFloat=new EigsFloat();
		static GaussJordanDouble statGaussJordanDouble=new GaussJordanDouble();
		static GaussJordanFloat statGaussJordanFloat=new GaussJordanFloat();
		static InverseDouble statInverseDouble=new InverseDouble();
		static Pythag statPythag=new Pythag();
		static SvdDouble statSvdDouble=new SvdDouble();
		static SvdFloat statSvdFloat=new SvdFloat();
		static Swap statSwap=new Swap();
		public delegate void PulseFunc();
		public delegate double DelChisq(double[] a);
		//	Hold static memory
		static Nr(){}
		/// <summary>
		/// Cholesky factorization.
		/// </summary>
		public class CholFloat
		{
			private static int i,j,k,n;
			private static float sum;
			public CholFloat() {}

			/// <summary>
			/// Computes the Cholesky factor of a symmetric, positive semidefinite matrix U.  If the matrix is not positive semidefinite, "false" is returned
			/// and "output" is meaningless.  However, if the matrix is positive semidefinite, then true is returned and "output" contains the Cholesky factor.
			/// </summary>
			/// <param name="U">input, The matrix to be factored.  Only the upper triangle is used (including the diagonal).  Lower triangle is ignored.</param>
			/// <param name="output">output, The Cholesky factor as a lower-triangular matrix.  Elements above the diagonal are meaningless.</param>
			/// <returns>False if the matrix was not positive semidefinite, meaning that a Cholesky factor could not be computed.  True otherwise.</returns>
			public static bool Run(float[,] U, float[,] output)
			{
				n = U.GetLength(0);
				for (i=0;i<n;i++)
				{
					for (j=i;j<n;j++)
					{
						sum=U[i,j];
						for (k=i-1;k>=0;k--)
							sum -= U[i,k]*U[j,k];
						if (i == j)
						{
							if (sum <= 0.0)
								return false;
							output[i,i]=(float)Math.Sqrt((double)sum);
						}
						else
							output[j,i]=sum/output[i,i];
					}
				}
				return true;
			}
		}
		/// <summary>
		/// Matrix inversion.
		/// </summary>
		public class InverseDouble
		{
			private static double d;
			private static double[,] a;
			private static int[] indx, b;
			private static double[] col, vv;
			private static double TINY = 1.0e-20;
			public InverseDouble() {}
			/// <summary>
			/// Computes the inverse of a matrix using LU decomposition and backsubstitution.
			/// </summary>
			/// <param name="A">The matrix to be evaluated for its inverse.</param>
			/// <param name="output">The output, which is the inverse of A.  A is not altered.</param>
			/// <returns>Returns true if successful.  False means A is likely to be singular.</returns>
			public static bool Run(double[,] A, double[,] output)
			{
				if(a==null || a.GetLength(0)<A.GetLength(0))
					a = new double[A.GetLength(0), A.GetLength(0)];
				int i,j;
				int len = A.GetLength(0);
				for(i=0; i<len; i++)
					for(j=0; j<len; j++)
						a[i,j]=A[i,j];

				if(indx==null || indx.Length < len)
					indx = new int[len];

				if(b==null || b.Length < len)
					b = new int[len];

				d = 1.0;

				if(col==null || col.Length < len)
					col = new double[len];
				
				if(!Ludcmp(a, len, indx, ref d))
					return false;
				else
				{
					for(j=0; j < len; j++)
					{
						for(i=0; i < len; i++)
							col[i] = 0.0;
						col[j] = 1.0;
						Lubksb(a, len, indx, col);
						for(i=0; i < len; i++)
							output[i,j] = col[i];
					}
				}
				return true;
			}
			/// <summary>
			/// Taken from Numerical Recipes.  Calculates a LU Decomposition of a column permutation of a.
			/// Outputs via reference to a and indx.  These values will be changed after using this function.
			/// </summary>
			/// <param name="a">The input matrix for LU Decomp.  Outputs the LU Decomp of a column 
			/// permutation of a.</param>
			/// <param name="n">The length of a side of the square matrix, 'a'.</param>
			/// <param name="indx">This is an output vector that records the row permutations created 
			/// by partial pivoting.</param>
			/// <param name="d">This is an output value that equals + or - 1.0 depending on whether or 
			/// no the number of
			/// row interchanges was even or odd respectively. </param>
			/// <returns>False if a is singular, true otherwise.</returns>
			public static bool Ludcmp(double[,] a, int n, int[] indx, ref double d)
			{
				int i,j,k;
				double big,dum,sum,temp;
				int imax = 0;
				
				if(vv==null || vv.Length<n)
					vv = new double[n];

				d=1.0;
				for (i=0;i<n;i++) 
				{
					big=0.0;
					for (j=0;j<n;j++)
						if((temp=Math.Abs(a[i,j])) > big) big=temp;
					if (big == 0.0)
						return false;
					vv[i]=1.0/big;
				}
				for (j=0;j<n;j++) 
				{
					for (i=0;i<j;i++) 
					{
						sum=a[i,j];
						for (k=0;k<i;k++)
							sum -= a[i,k]*a[k,j];
						a[i,j]=sum;
					}
					big=0.0;
					for (i=j;i<n;i++) 
					{
						sum=a[i,j];
						for (k=0;k<j;k++)
							sum -= a[i,k]*a[k,j];
						a[i,j]=sum;
						if ( (dum=vv[i]*Math.Abs(sum)) >= big) 
						{
							big=dum;
							imax=i;
						}
					}
					if (j != imax) 
					{
						for (k=0;k<n;k++) 
						{
							dum=a[imax,k];
							a[imax,k]=a[j,k];
							a[j,k]=dum;
						}
						d = -d;
						vv[imax]=vv[j];
					}
					indx[j]=imax;
					if (a[j,j] == 0.0) a[j,j]=TINY;
					if (j != n-1) 
					{
						dum=1.0/(a[j,j]);
						for (i=j+1;i<n;i++)
							a[i,j] *= dum;
					}
				}
				return true;
			}

			/// <summary>
			/// Taken from Numerical Recipes in C.
			/// Uses LU Back-Substitution that solves the set of n linear equations AX=B.
			/// </summary>
			/// <param name="a">An LU Decomposition of "A", determined by the function ludcmp.</param>
			/// <param name="n">The length of a side of the square matrix 'a'.</param>
			/// <param name="indx">The index into the row permutations of the LU Decomposition.</param>
			/// <param name="b">The variable "B", also used for output.</param>
			/// <returns>"X", the solution to AX=B</returns>
			public static void Lubksb(double[,] a, int n, int[] indx, double[] b)
			{
				int i,ii=0,ip,j;
				double sum;

				for (i=0;i<n;i++) 
				{
					ip=indx[i];
					sum=b[ip];
					b[ip]=b[i];
					if (ii != 0)
						for (j=ii-1;j<i;j++)
							sum -= a[i,j]*b[j];
					else if (sum != 0.0)
						ii=i+1;
					b[i]=sum;
				}
				for (i=n-1;i>=0;i--) 
				{
					sum=b[i];
					for (j=i+1;j<n;j++)
						sum -= a[i,j]*b[j];
					b[i]=sum/a[i,i];
				}
			}
		}
		/// <summary>
		/// Eigensystem determination.
		/// </summary>
		public class EigsDouble
		{
			private static double scale,hh,h,g,f,s,r,p,dd,c,b;
			private static int l,k,j,i, m, iter;
			private static double[] _e;
			/// <summary>
			/// Computes the eigensystem of mat[0..n-1,0..n-1] without changing it.  Outputs eigenvalues as 'd'
			/// and corresponding eigenvectors as columns of 'Q'.  Information is sorted by eigenvalue.
			/// </summary>
			/// <param name="mat"></param>
			/// <param name="n"></param>
			/// <param name="d"></param>
			/// <param name="Q"></param>
			public static void RunRealSymmetric(double[,] mat,  int n, double[] d, double[,] Q)
			{
				//	Copy the values of mat to Q.
				for(i=0; i<n; i++)
					for(j=0; j<n; j++)
						Q[i,j] = mat[i,j];
				//	Run the algorithm.
				RunRealSymmetric(Q, n, d);
			}
			/// <summary>
			/// Computes the eigensystem of inMatOutQ[0..n-1,0..n-1].  Outputs eigenvalues as 'd'
			/// and corresponding eigenvectors as columns of inMatOutQ[0..n-1,0..n-1].  Information is sorted by eigenvalue.
			/// </summary>
			/// <param name="inMatOutQ"></param>
			/// <param name="n"></param>
			/// <param name="d"></param>
			public static void RunRealSymmetric(double[,] inMatOutQ, int n, double[] d)
			{
				//	Static resources
				if(_e==null || _e.Length<n)
					_e = new double[n];
				//	Convert to symmetric tridiagonal form.
				Tred2(inMatOutQ, n, d, _e);
				//	QL algorithm with implicit shifts.
				Tqli(n, d, _e, inMatOutQ);
				//	Sort the results.
				Eigsrt(d, inMatOutQ, n);
			}
			/// <summary>
			/// Sorts the eigenvalues and corresponding eigenvector columns.
			/// </summary>
			/// <param name="d">Eigenvalues.</param>
			/// <param name="v">Columns of eigenvectors.</param>
			/// <param name="n">Number of eigenvalues and eigenvectors.</param>
			public static void Eigsrt(double[] d, double[,] v, int n)
			{
				for (i=0;i<n-1;i++)
				{
					p=d[k=i];
					for (j=i;j<n;j++)
						if (d[j] >= p) p=d[k=j];
					if (k != i)
					{
						d[k]=d[i];
						d[i]=p;
						for (j=0;j<n;j++)
						{
							p=v[j,i];
							v[j,i]=v[j,k];
							v[j,k]=p;
						}
					}
				}
			}

			/// <summary>
			/// QL algorithm with implicit shifts, to determine the eigenvalues and eigenvectors of a real, symmetric,
			/// tridiagonal matrix, or of a real, symmetric matrix previously reduced by tred2.
			/// </summary>
			/// <param name="n">int  The size of the vector of eigenvalues.</param>
			/// <param name="d">In/Out, double[], On input, d[0..n-1] gives the diagonal elements of the tridiagonal matrix.
			/// On, output, it contains the eigenvalues.</param>
			/// <param name="e">In/Out, The sub-diagonal elements of the tridiagonal matrix, destroyed on output.</param>
			/// <param name="z">In/Out, On input, this contains the matrix that is output by  the Tred2 function (argument 'a'),
			/// also known as the Q-matrix.  On output, this contains the eigenvectors associated with the 'd' vector.</param>
			public static void Tqli(int n, double[] d, double[] e, double[,] z)
			{
				for (i=1;i<n;i++)
					e[i-1]=e[i];
				e[n-1]=0.0;
				for (l=0;l<n;l++)
				{
					iter=0;
					do
					{
						for (m=l;m<n-1;m++)
						{
							dd=Math.Abs(d[m])+Math.Abs(d[m+1]);
							if (Math.Abs(e[m])+dd == dd) break;
						}
						if (m != l)
						{
							if (iter++ == 30)
								throw new Exception("Too many iterations in Tqli while the eigensystem is being determined.");
							g=(d[l+1]-d[l])/(2.0*e[l]);
							r=Pythag.Run(g,1.0);
							g=d[m]-d[l]+e[l]/(g+Sign(r,g));
							s=c=1.0;
							p=0.0;
							for (i=m-1;i>=l;i--)
							{
								f=s*e[i];
								b=c*e[i];
								e[i+1]=(r=Pythag.Run(f,g));
								if (r == 0.0)
								{
									d[i+1] -= p;
									e[m]=0.0;
									break;
								}
								s=f/r;
								c=g/r;
								g=d[i+1]-p;
								r=(d[i]-g)*s+2.0*c*b;
								d[i+1]=g+(p=s*r);
								g=c*r-b;
								// Next loop can be omitted if eigenvectors not wanted
								for (k=0;k<n;k++)
								{
									f=z[k,i+1];
									z[k,i+1]=s*z[k,i]+c*f;
									z[k,i]=c*z[k,i]-s*f;
								}
							}
							if (r == 0.0 && i >= l)
								continue;
							d[l] -= p;
							e[l]=g;
							e[m]=0.0;
						}
					}
					while (m != l);
				}
			}
			/// <summary>
			/// Householder reduction of a real, symmetric matrix a[0..n-1,0..n-1].  On output, 'a' is replaced by the orthogonal
			/// matrix Q effecting the transformation.
			/// </summary>
			/// <param name="a">In/Out, Input as the matrix to be reduced, changed into an orthogonal matrix Q effecting the transformation.</param>
			/// <param name="n">The size of a side of the square matrix 'a'</param>
			/// <param name="d">Out  Diagonal elements of the tridiagonal matrix.</param>
			/// <param name="e">Out  Off-diagonal elements of the tridiagonal matrix.</param>
			public static void Tred2(double[,] a, int n, double[] d, double[] e)
			{
				for (i=n-1;i>0;i--)
				{
					l=i-1;
					h=scale=0.0;
					if (l > 0)
					{
						for (k=0;k<l+1;k++)
							scale += Math.Abs(a[i,k]);
						if (scale == 0.0)
							e[i]=a[i,l];
						else
						{
							for (k=0;k<l+1;k++)
							{
								a[i,k] /= scale;
								h += a[i,k]*a[i,k];
							}
							f=a[i,l];
							g=(f >= 0.0 ? -Math.Sqrt(h) : Math.Sqrt(h));
							e[i]=scale*g;
							h -= f*g;
							a[i,l]=f-g;
							f=0.0;
							for (j=0;j<l+1;j++)
							{
								// Next statement can be omitted if eigenvectors not wanted
								a[j,i]=a[i,j]/h;
								g=0.0;
								for (k=0;k<j+1;k++)
									g += a[j,k]*a[i,k];
								for (k=j+1;k<l+1;k++)
									g += a[k,j]*a[i,k];
								e[j]=g/h;
								f += e[j]*a[i,j];
							}
							hh=f/(h+h);
							for (j=0;j<l+1;j++)
							{
								f=a[i,j];
								e[j]=g=e[j]-hh*f;
								for (k=0;k<j+1;k++)
									a[j,k] -= (f*e[k]+g*a[i,k]);
							}
						}
					}
					else
						e[i]=a[i,l];
					d[i]=h;
				}
				// Next statement can be omitted if eigenvectors not wanted
				d[0]=0.0;
				e[0]=0.0;
				// Contents of this loop can be omitted if eigenvectors not
				//	wanted except for statement d[i]=a[i,i];
				for (i=0;i<n;i++)
				{
					l=i;
					if (d[i] != 0.0)
					{
						for (j=0;j<l;j++)
						{
							g=0.0;
							for (k=0;k<l;k++)
								g += a[i,k]*a[k,j];
							for (k=0;k<l;k++)
								a[k,j] -= g*a[k,i];
						}
					}
					d[i]=a[i,i];
					a[i,i]=1.0;
					for(j=0;j<l;j++)
						a[j,i]=a[i,j]=0.0;
				}
			}
		}
		/// <summary>
		/// Eigensystem determination.
		/// </summary>
		public class EigsFloat
		{
			private static float scale,hh,h,g,f,s,r,p,dd,c,b;
			private static int l,k,j,i, m, iter;
			private static float[] _e;
			/// <summary>
			/// Computes the eigensystem of mat[0..n-1,0..n-1] without changing it.  Outputs eigenvalues as 'd'
			/// and corresponding eigenvectors as columns of 'Q'.  Information is sorted by eigenvalue.
			/// </summary>
			/// <param name="mat"></param>
			/// <param name="n"></param>
			/// <param name="d"></param>
			/// <param name="Q"></param>
			public static void RunRealSymmetric(float[,] mat,  int n, float[] d, float[,] Q)
			{
				//	Copy the values of mat to Q.
				for(i=0; i<n; i++)
					for(j=0; j<n; j++)
						Q[i,j] = mat[i,j];
				//	Run the algorithm.
				RunRealSymmetric(Q, n, d);
			}
			/// <summary>
			/// Computes the eigensystem of inMatOutQ[0..n-1,0..n-1].  Outputs eigenvalues as 'd'
			/// and corresponding eigenvectors as columns of inMatOutQ[0..n-1,0..n-1].  Information is sorted by eigenvalue.
			/// </summary>
			/// <param name="inMatOutQ"></param>
			/// <param name="n"></param>
			/// <param name="d"></param>
			public static void RunRealSymmetric(float[,] inMatOutQ, int n, float[] d)
			{
				//	Static resources
				if(_e==null || _e.Length<n)
					_e = new float[n];
				//	Convert to symmetric tridiagonal form.
				Tred2(inMatOutQ, n, d, _e);
				//	QL algorithm with implicit shifts.
				Tqli(n, d, _e, inMatOutQ);
				//	Sort the results.
				Eigsrt(d, inMatOutQ, n);
			}
			/// <summary>
			/// Sorts the eigenvalues and corresponding eigenvector columns.
			/// </summary>
			/// <param name="d">Eigenvalues.</param>
			/// <param name="v">Columns of eigenvectors.</param>
			/// <param name="n">Number of eigenvalues and eigenvectors.</param>
			public static void Eigsrt(float[] d, float[,] v, int n)
			{
				for (i=0;i<n-1;i++)
				{
					p=d[k=i];
					for (j=i;j<n;j++)
						if (d[j] >= p) p=d[k=j];
					if (k != i)
					{
						d[k]=d[i];
						d[i]=p;
						for (j=0;j<n;j++)
						{
							p=v[j,i];
							v[j,i]=v[j,k];
							v[j,k]=p;
						}
					}
				}
			}

			/// <summary>
			/// QL algorithm with implicit shifts, to determine the eigenvalues and eigenvectors of a real, symmetric,
			/// tridiagonal matrix, or of a real, symmetric matrix previously reduced by tred2.
			/// </summary>
			/// <param name="n">int  The size of the vector of eigenvalues.</param>
			/// <param name="d">In/Out, float[], On input, d[0..n-1] gives the diagonal elements of the tridiagonal matrix.
			/// On, output, it contains the eigenvalues.</param>
			/// <param name="e">In/Out, The sub-diagonal elements of the tridiagonal matrix, destroyed on output.</param>
			/// <param name="z">In/Out, On input, this contains the matrix that is output by  the Tred2 function (argument 'a'),
			/// also known as the Q-matrix.  On output, this contains the eigenvectors associated with the 'd' vector.</param>
			public static void Tqli(int n, float[] d, float[] e, float[,] z)
			{
				for (i=1;i<n;i++)
					e[i-1]=e[i];
				e[n-1]=0.0f;
				for (l=0;l<n;l++)
				{
					iter=0;
					do
					{
						for (m=l;m<n-1;m++)
						{
							dd=Math.Abs(d[m])+Math.Abs(d[m+1]);
							if (Math.Abs(e[m])+dd == dd) break;
						}
						if (m != l)
						{
							if (iter++ == 30)
								throw new Exception("Too many iterations in Tqli while the eigensystem is being determined.");
							g=(d[l+1]-d[l])/(2.0f*e[l]);
							r=Pythag.Run(g,1.0f);
							g=d[m]-d[l]+e[l]/(g+Sign(r,g));
							s=c=1.0f;
							p=0.0f;
							for (i=m-1;i>=l;i--)
							{
								f=s*e[i];
								b=c*e[i];
								e[i+1]=(r=Pythag.Run(f,g));
								if (r == 0.0f)
								{
									d[i+1] -= p;
									e[m]=0.0f;
									break;
								}
								s=f/r;
								c=g/r;
								g=d[i+1]-p;
								r=(d[i]-g)*s+2.0f*c*b;
								d[i+1]=g+(p=s*r);
								g=c*r-b;
								// Next loop can be omitted if eigenvectors not wanted
								for (k=0;k<n;k++)
								{
									f=z[k,i+1];
									z[k,i+1]=s*z[k,i]+c*f;
									z[k,i]=c*z[k,i]-s*f;
								}
							}
							if (r == 0.0 && i >= l)
								continue;
							d[l] -= p;
							e[l]=g;
							e[m]=0.0f;
						}
					}
					while (m != l);
				}
			}
			/// <summary>
			/// Householder reduction of a real, symmetric matrix a[0..n-1,0..n-1].  On output, 'a' is replaced by the orthogonal
			/// matrix Q effecting the transformation.
			/// </summary>
			/// <param name="a">In/Out, Input as the matrix to be reduced, changed into an orthogonal matrix Q effecting the transformation.</param>
			/// <param name="n">The size of a side of the square matrix 'a'</param>
			/// <param name="d">Out  Diagonal elements of the tridiagonal matrix.</param>
			/// <param name="e">Out  Off-diagonal elements of the tridiagonal matrix.</param>
			public static void Tred2(float[,] a, int n, float[] d, float[] e)
			{
				for (i=n-1;i>0;i--)
				{
					l=i-1;
					h=scale=0.0f;
					if (l > 0)
					{
						for (k=0;k<l+1;k++)
							scale += Math.Abs(a[i,k]);
						if (scale == 0.0f)
							e[i]=a[i,l];
						else
						{
							for (k=0;k<l+1;k++)
							{
								a[i,k] /= scale;
								h += a[i,k]*a[i,k];
							}
							f=a[i,l];
							g=(f >= 0.0f ? -(float)Math.Sqrt((double)h) : (float)Math.Sqrt((double)h));
							e[i]=scale*g;
							h -= f*g;
							a[i,l]=f-g;
							f=0.0f;
							for (j=0;j<l+1;j++)
							{
								// Next statement can be omitted if eigenvectors not wanted
								a[j,i]=a[i,j]/h;
								g=0.0f;
								for (k=0;k<j+1;k++)
									g += a[j,k]*a[i,k];
								for (k=j+1;k<l+1;k++)
									g += a[k,j]*a[i,k];
								e[j]=g/h;
								f += e[j]*a[i,j];
							}
							hh=f/(h+h);
							for (j=0;j<l+1;j++)
							{
								f=a[i,j];
								e[j]=g=e[j]-hh*f;
								for (k=0;k<j+1;k++)
									a[j,k] -= (f*e[k]+g*a[i,k]);
							}
						}
					}
					else
						e[i]=a[i,l];
					d[i]=h;
				}
				// Next statement can be omitted if eigenvectors not wanted
				d[0]=0.0f;
				e[0]=0.0f;
				// Contents of this loop can be omitted if eigenvectors not
				//	wanted except for statement d[i]=a[i,i];
				for (i=0;i<n;i++)
				{
					l=i;
					if (d[i] != 0.0)
					{
						for (j=0;j<l;j++)
						{
							g=0.0f;
							for (k=0;k<l;k++)
								g += a[i,k]*a[k,j];
							for (k=0;k<l;k++)
								a[k,j] -= g*a[k,i];
						}
					}
					d[i]=a[i,i];
					a[i,i]=1.0f;
					for(j=0;j<l;j++)
						a[j,i]=a[i,j]=0.0f;
				}
			}
		}
		/// <summary>
		/// Facilitates Gauss Jordan elimination.
		/// </summary>
		public class GaussJordanDouble
		{
			static int[] indxc,indxr,ipiv;
			static double big,dum,pivinv;
			/// <summary>
			/// Linear equation solution by Gauss-Jordan elimination.
			/// </summary>
			/// <param name="a">double[1..n,1..n] Input matrix.  On output, this is set to the
			/// inverse of the input.</param>
			/// <param name="n">dimension of the square input matrix</param>
			/// <param name="b">double[1..n,1..m] The right hand side vectors.  On output, this is
			/// replaced by the corresponding set of solution vectors.</param>
			/// <param name="m">number of right hand side vectors in b</param>
			public static void Run(double[,] a, int n, double[,] b, int m)
			{
				if(indxc==null || indxc.Length<n)
				{
					indxc=new int[n];
					indxr=new int[n];
					ipiv=new int[n];
				}

				int i,icol=0,irow=0,j,k,l,ll;
				
				for(j=0;j<n;j++)
					ipiv[j]=0;
				for(i=0;i<n;i++) 
				{
					big=0.0;
					for (j=0;j<n;j++)
						if (ipiv[j] != 1)
							for (k=0;k<n;k++) 
							{
								if (ipiv[k] == 0) 
								{
									if (Math.Abs(a[j,k]) >= big) 
									{
										big=Math.Abs(a[j,k]);
										irow=j;
										icol=k;
									}
								} 
							}
					ipiv[icol]++;
					if (irow != icol) 
					{
						for (l=0;l<n;l++)
							Swap.Run(ref a[irow,l],ref a[icol,l]);
						for (l=0;l<m;l++)
							Swap.Run(ref b[irow,l],ref b[icol,l]);
					}
					indxr[i]=irow;
					indxc[i]=icol;
					if (a[icol,icol] == 0.0)
						throw new System.Exception("MathClasses.Matrix.GaussJ:  Singular Matrix.");;
					pivinv=1.0/a[icol,icol];
					a[icol,icol]=1.0;
					for(l=0;l<n;l++)
						a[icol,l] *= pivinv;
					for(l=0;l<m;l++)
						b[icol,l] *= pivinv;
					for(ll=0;ll<n;ll++)
						if (ll != icol) 
						{
							dum=a[ll,icol];
							a[ll,icol]=0.0;
							for (l=0;l<n;l++)
								a[ll,l] -= a[icol,l]*dum;
							for (l=0;l<m;l++)
								b[ll,l] -= b[icol,l]*dum;
						}
				}
				for(l=n-1;l>=0;l--)
					if(indxr[l] != indxc[l])
						for(k=0;k<n;k++)
							Swap.Run(ref a[k,indxr[l]],ref a[k,indxc[l]]);
			}
		}
		/// <summary>
		/// Facilitates Gauss Jordan elimination.
		/// </summary>
		public class GaussJordanFloat
		{
			static int[] indxc,indxr,ipiv;
			static float big,dum,pivinv;
			/// <summary>
			/// Linear equation solution by Gauss-Jordan elimination.
			/// </summary>
			/// <param name="a">float[1..n,1..n] Input matrix.  On output, this is set to the
			/// inverse of the input.</param>
			/// <param name="n">dimension of the square input matrix</param>
			/// <param name="b">float[1..n,1..m] The right hand side vectors.  On output, this is
			/// replaced by the corresponding set of solution vectors.</param>
			/// <param name="m">number of right hand side vectors in b</param>
			public static void Run(float[,] a, int n, float[,] b, int m)
			{
				if(indxc==null || indxc.Length<n)
				{
					indxc=new int[n];
					indxr=new int[n];
					ipiv=new int[n];
				}

				int i,icol=0,irow=0,j,k,l,ll;
				
				for(j=0;j<n;j++)
					ipiv[j]=0;
				for(i=0;i<n;i++) 
				{
					big=0.0f;
					for (j=0;j<n;j++)
						if (ipiv[j] != 1)
							for (k=0;k<n;k++) 
							{
								if (ipiv[k] == 0) 
								{
									if (Math.Abs(a[j,k]) >= big) 
									{
										big=Math.Abs(a[j,k]);
										irow=j;
										icol=k;
									}
								} 
							}
					ipiv[icol]++;
					if (irow != icol) 
					{
						for (l=0;l<n;l++)
							Swap.Run(ref a[irow,l],ref a[icol,l]);
						for (l=0;l<m;l++)
							Swap.Run(ref b[irow,l],ref b[icol,l]);
					}
					indxr[i]=irow;
					indxc[i]=icol;
					if (a[icol,icol] == 0.0)
						throw new System.Exception("MathClasses.Matrix.GaussJ:  Singular Matrix.");;
					pivinv=1.0f/a[icol,icol];
					a[icol,icol]=1.0f;
					for(l=0;l<n;l++)
						a[icol,l] *= pivinv;
					for(l=0;l<m;l++)
						b[icol,l] *= pivinv;
					for(ll=0;ll<n;ll++)
						if (ll != icol) 
						{
							dum=a[ll,icol];
							a[ll,icol]=0.0f;
							for (l=0;l<n;l++)
								a[ll,l] -= a[icol,l]*dum;
							for (l=0;l<m;l++)
								b[ll,l] -= b[icol,l]*dum;
						}
				}
				for(l=n-1;l>=0;l--)
					if(indxr[l] != indxc[l])
						for(k=0;k<n;k++)
							Swap.Run(ref a[k,indxr[l]],ref a[k,indxc[l]]);
			}
		}
		/// <summary>
		/// Expand in storage the covariance matrix "covar" to take into account the fixed parameters.
		/// </summary>
		/// <param name="covar">double[,] Covariance matrix.</param>
		/// <param name="ma"></param>
		/// <param name="ia"></param>
		/// <param name="mfit"></param>
		public static void CovSrt(double[,] covar, int ma, bool[] ia, int mfit)
		{
			int i,j,k;
			for(i=mfit;i<ma;i++)
				for(j=0;j<i+1;j++)
					covar[i,j]=covar[j,i]=0.0;
			k=mfit-1;
			for (j=ma-1;j>=0;j--)
			{
				if (ia[j]) 
				{
					for (i=0;i<ma;i++)
						Swap.Run(ref covar[i,k],ref covar[i,j]);
					for (i=0;i<ma;i++)
						Swap.Run(ref covar[k,i],ref covar[j,i]);
					k--;
				}
			}
		}
		/// <summary>
		/// Expand in storage the covariance matrix "covar" to take into account the fixed parameters.
		/// </summary>
		/// <param name="covar">double[,] Covariance matrix.</param>
		/// <param name="ma"></param>
		/// <param name="ia"></param>
		/// <param name="mfit"></param>
		public static void CovSrt(float[,] covar, int ma, bool[] ia, int mfit)
		{
			int i,j,k;
			for(i=mfit;i<ma;i++)
				for(j=0;j<i+1;j++)
					covar[i,j]=covar[j,i]=0.0f;
			k=mfit-1;
			for (j=ma-1;j>=0;j--)
			{
				if (ia[j]) 
				{
					for (i=0;i<ma;i++)
						Swap.Run(ref covar[i,k],ref covar[i,j]);
					for (i=0;i<ma;i++)
						Swap.Run(ref covar[k,i],ref covar[j,i]);
					k--;
				}
			}
		}
		/// <summary>
		/// Minimization algorithms taken from Numerical Recipes in C++.  Altered for C#.
		/// Put into a class for optimal memory management on repeated usage.
		/// </summary>
		public class LevMarDouble
		{
			/// <summary>
			/// For conditioning SVD in MrqMin()
			/// </summary>
			private double svd_wmax;
			/// <summary>
			/// Used as working space during iterations.  On finish, gives the covariance matrix of noise for the free parameters.
			/// </summary>
			public double[,] Covar;
			/// <summary>
			/// Used as working space during iterations.  On finish, gives the curvature matrix of free parameters at the solution point.
			/// </summary>
			public double[,] Alpha;
			public double[] Beta;
			/// <summary>
			/// Used as working space during iterations in MrqCof.
			/// </summary>
			public double[] Dyda;
			private double chiTry, ymod, wt, sig2i, dy, chiImprovement;
			private double[] atry, da, onedaTemp, onedaTempOut, wVec;
			private double[,] covarTemp, Vmatrix, oneda;
			/// <summary>
			/// Allocates the internal memory necessary for execution (if necessary).  Called locally in Run().
			/// </summary>
			/// <param name="nParams"></param>
			public void InternalMalloc(int nParams)
			{
				if(atry==null  ||  nParams > atry.Length)
				{
					atry=new double[nParams];
					Beta=new double[nParams];
					da=new double[nParams];

					Alpha = new double[nParams, nParams];
					Covar = new double[nParams, nParams];
					Dyda = new double[nParams];

					//	These were supposed to be Nfit length instead of nParams.  nParams is always longer or equal to Nfit.
					//	I changed this so that static memory was always determined by nParams, and Nfit has nothing to do with it.
					oneda = new double[nParams,1];
					onedaTemp = new double[nParams];
					onedaTempOut = new double[nParams];
					Vmatrix = new double[nParams,nParams];
					wVec = new double[nParams];
					covarTemp = new double[nParams,nParams];
				}
			}
			/// <summary>
			/// The X values.
			/// </summary>
			public double[] X;
			/// <summary>
			/// The Y values.
			/// </summary>
			public double[] Y;
			/// <summary>
			/// The standard deviation associated with each Y-->f(X) observation in the data.
			/// </summary>
			public double[] Sd;
			/// <summary>
			/// The length of the data set, or the number of significant elements in X, Y, and Sd.
			/// </summary>
			public int Ndata;
			/// <summary>
			/// An improvement in Chisq of this small a value has no meaning.  Thus, iteration stops when improvements are this small.
			/// </summary>
			public double ChisqErrTolerance;
			/// <summary>
			/// This delegates a function that calculates the predicted value, y-->f(x), for a given set of parameters, a.
			/// This function also gives the derivative of y with respect to changes in a.
			/// </summary>
			public delegate double MrqFunc(double x, double[] a, double[] dyda);
			/// <summary>
			/// The model which takes, as arguments, a set of free parameters in order to generate predictions of the model.
			/// </summary>
			public MrqFunc Func;
			/// <summary>
			/// The boolean vector that determines whether a particular free parameter (represented by A) needs to be optimized.
			/// If this is false, the parameter remains constant during the optimization procedure.
			/// </summary>
			public bool[] Op;
			/// <summary>
			/// The free parameters for the model.
			/// </summary>
			public double[] A;
			/// <summary>
			/// A Hessian-Jacobian tradeoff factor.
			/// </summary>
			public double Lamda;
			/// <summary>
			/// The value of chi-squared for the set of free parameters (A).
			/// </summary>
			public double Chisq;
			//	The number of true values in the Op vector.  Updated in the constructor.
			public int Nfit;
			/// <summary>
			/// Counts the number of iterations through MrqMin();
			/// </summary>
			public int Iterations;
			/// <summary>
			/// Run the minimization routine.
			/// </summary>
			/// <param name="Ndata"></param>
			/// <param name="X"></param>
			/// <param name="Y"></param>
			/// <param name="Sd"></param>
			/// <param name="Func"></param>
			/// <param name="A"></param>
			/// <param name="Op"></param>
			/// <param name="ChisqErrTolerance"></param>
			public void Run(int Ndata, double[] X, double[] Y, double[] Sd, MrqFunc Func, double[] A, bool[] Op, double ChisqErrTolerance)
			{
				int i;

				//	Set the member variables accordingly.
				this.X=X;
				this.Y=Y;
				this.Sd=Sd;
				this.Func = Func;
				this.A = A;
				this.Op = Op;
				this.ChisqErrTolerance = ChisqErrTolerance;
				this.Ndata = Ndata;

				//	Insure that static memory is properly allocated.
				InternalMalloc(A.Length);

				//	Compute Nfit;
				Nfit = 0;
				for(i=0; i<Op.Length; i++)
					if(Op[i])
						Nfit++;

				//	Set Lamda to 0.001 as an initial value.  Conduct the first fit.				
				Lamda=0.001;
				//	First pass, MrqCof
				Iterations=0;
				Chisq = MrqCof(A, Alpha, Beta);
				//	Copy the values to atry, the vector of free parameters to try.
				A.CopyTo(atry, 0);
				//	Go
				MrqMin();
			}
			/// <summary>
			/// Initiates a Levenberg-Marquardt fit of the data provided.  On output, the parameters A are changed to their optimal values that
			/// have minimized sum-squared error.  The member variable 
			/// </summary>
			/// <param name="Ndata">int  The length of the data set, or the number of significant elements in X, Y, and Sd.</param>
			/// <param name="X">double[]  Input only.  The independent variables.  This is not changed.</param>
			/// <param name="Y">double[]  Input only.  The data that is dependent on the independent variables.  This is not changed.</param>
			/// <param name="Sd">double[]  Input only.  The standard deviation of Y.  This is not changed.</param>
			/// <param name="Func">MrqFunc  Input only.  The model that generates a predicted value of Y based on a value of X and a set of free parameters A.</param>
			/// <param name="A">double[]  Input and output.  First guess at free parameters.  This becomes the member variable A.  These values are altered.</param>
			/// <param name="Op">double[]  Input only.  This becomes the member variable Op.  This is not altered.</param>
			/// <param name="ChisqErrTolerance">double greater than zero.
			/// An improvement in Chisq of this small a value has no meaning.  Thus, iteration stops when improvements are this small.</param>
			public LevMarDouble(int Ndata, double[] X, double[] Y, double[] Sd, MrqFunc Func, double[] A, bool[] Op, double ChisqErrTolerance)
			{
				Run(Ndata, X, Y, Sd, Func, A, Op, ChisqErrTolerance);
			}
			
			/// <summary>
			/// Used by MrqMin to evaluate the linearized fitting matrix "Alpha", vector "Beta".
			/// </summary>
			/// <returns>The value of chi-squared for this set of free parameters.</returns>
			public double MrqCof(double[] a, double[,] alpha, double[] beta)
			{
				int i,j,k,l,m;
				double chisq=0.0;
				
				for (j=0;j<Nfit;j++)
				{
					for (k=0;k<=j;k++)
						alpha[j,k]=0.0;
					beta[j]=0.0;
				}
				for (i=0;i<Ndata;i++)
				{
					ymod=Func(X[i],a,Dyda);
					sig2i=1.0/(Sd[i]*Sd[i]);
					dy=Y[i]-ymod;
					for (j=0,l=0;l<a.Length;l++)
					{
						if (Op[l])
						{
							wt=Dyda[l]*sig2i;
							for(k=0,m=0;m<l+1;m++)
								if(Op[m])
									alpha[j,k++] += wt*Dyda[m];
							beta[j++] += dy*wt;
						}
					}
					chisq += dy*dy*sig2i;
				}
				for(j=1;j<Nfit;j++)
					for(k=0;k<j;k++)
						alpha[k,j]=alpha[j,k];
				return chisq;
			}
			/// <summary>
			/// Drives the minimization routine.
			/// </summary>
			public void MrqMin()
			{
				Iterations++;
				int j,k,l;

				for(j=0;j<Nfit;j++)
				{
					for(k=0;k<Nfit;k++)
					{
						covarTemp[j,k]=Covar[j,k]=Alpha[j,k];
					}
					covarTemp[j,j]=Covar[j,j]=Alpha[j,j]*(1.0+Lamda);
					onedaTemp[j]=oneda[j,0]=Beta[j];
				}

				//	The following try/catch loop is an edition by Dave Ing.  The original algorithm is just "GaussJordan.Run(Covar,Nfit,oneda,1);"
				//	Use GaussJordan to invert the covariance matrix as fast as possible in addition to computing the optimal step.
				try
				{
					GaussJordanDouble.Run(Covar,Nfit,oneda,1);
				}
					//	Singularity, use SVD instead.
				catch
				{
					SvdDouble.Decomp(covarTemp, Nfit, Nfit, wVec, Vmatrix);
					svd_wmax=0.0;
					for(j=0; j<Nfit; j++)
						if(wVec[j]>svd_wmax)
							svd_wmax = wVec[j];
					svd_wmax/=50.0;
					for(j=0; j<Nfit; j++)
						if(wVec[j]<svd_wmax)
							wVec[j]=0.0;
					SvdDouble.Backsub(covarTemp, Nfit, Nfit, wVec, Vmatrix, onedaTemp, onedaTempOut);
					for(j=0; j<Nfit; j++)
						for(k=0; k<Nfit; k++)
							Covar[j,k]=covarTemp[j,k];
					for(j=0; j<Nfit; j++)
						oneda[j,0] = onedaTempOut[j];
				}
				
				for (j=0;j<Nfit;j++)
					da[j]=oneda[j,0];
				if (Lamda == 0.0)
				{
					CovSrt(Covar,A.Length,Op,Nfit);
					CovSrt(Alpha,A.Length,Op,Nfit);
					return;
				}
				for (j=0,l=0;l<A.Length;l++)
					if (Op[l])
						atry[l]=A[l]+da[j++];
				chiTry=MrqCof(atry, Covar, da);
				if (chiTry < Chisq)
				{
					Lamda *= 0.1;
					chiImprovement=Chisq-chiTry;
					Chisq=chiTry;
					for (j=0;j<Nfit;j++)
					{
						for (k=0;k<Nfit;k++)
							Alpha[j,k]=Covar[j,k];
						Beta[j]=da[j];
					}
					for (l=0;l<A.Length;l++) A[l]=atry[l];
					//	Stopping criteria:  Is the Chisq improvement less than the ErrTolerance parameter.
					if(chiImprovement < ChisqErrTolerance)
						Lamda=0.0;	//	stop algorithm on next run.
					MrqMin();
				}
				else
				{
					Lamda *= 10.0;
					if(Lamda>1e10)
						Lamda=0.0;
					//	Reiterate on failure unless Lamda is very large.
					MrqMin();
				}
			}
			/// <summary>
			/// Construct an empty object.
			/// </summary>
			public LevMarDouble() {}
		}
		public class SvdDouble
		{
			private static double anorm,c,f,g,h,s,scale,x,y,z;
			private static double[] rv1, tmp;
			private static int i,its,j,jj,k,l=0,nm=0;
			private static bool flag;
			/// <summary>
			/// Performs singular value back-substitution into x (A=U*diag(w)*V') after Svdcmp.
			/// </summary>
			/// <param name="u">Input matrix, the "a" argument after Svdcmp was run on it.</param>
			/// <param name="m">The number of rows in a.</param>
			/// <param name="n">The number of columns in a.</param>
			/// <param name="w">Output vector from Svdcmp.</param>
			/// <param name="v">Output vector from Svdcmp.</param>
			/// <param name="b">Right hand side input.</param>
			/// <param name="x">Solution vector output.</param>
			public static void Backsub(double[,] u, int m, int n, double[] w, double[,] v, double[] b, double[] x)
			{
				if(tmp==null || tmp.Length<n)
					tmp=new double[n];
				for (j=0;j<n;j++) 
				{
					s=0.0;
					if (w[j] != 0.0) 
					{
						for (i=0;i<m;i++) s += u[i,j]*b[i];
						s /= w[j];
					}
					tmp[j]=s;
				}
				for (j=0;j<n;j++) 
				{
					s=0.0;
					for (jj=0;jj<n;jj++) s += v[j,jj]*tmp[jj];
					x[j]=s;
				}
			}
			/// <summary>
			/// Performs the Singular Value Decomposition of "a" (matrix A), which is replaced by matrix U
			/// on output (for input into Svdbks).  w is the diagonal component of W.  v is V.  A=UWV'.
			/// </summary>
			/// <param name="a"></param>
			/// <param name="m">int  The number of rows in a.</param>
			/// <param name="n">int  The number of columns in a.</param>
			/// <param name="w"></param>
			/// <param name="v"></param>
			public static void Decomp(double[,] a, int m, int n, double[] w, double[,] v)
			{	
				if(rv1==null || rv1.Length<n)
					rv1=new double[n];
				g=scale=anorm=0.0;
				for (i=0;i<n;i++) 
				{
					l=i+2;
					rv1[i]=scale*g;
					g=s=scale=0.0;
					if (i < m) 
					{
						for (k=i;k<m;k++) scale += Math.Abs(a[k,i]);
						if (scale != 0.0) 
						{
							for (k=i;k<m;k++) 
							{
								a[k,i] /= scale;
								s += a[k,i]*a[k,i];
							}
							f=a[i,i];
							g = -Sign(Math.Sqrt(s),f);
							h=f*g-s;
							a[i,i]=f-g;
							for (j=l-1;j<n;j++) 
							{
								for (s=0.0,k=i;k<m;k++) s += a[k,i]*a[k,j];
								f=s/h;
								for (k=i;k<m;k++) a[k,j] += f*a[k,i];
							}
							for (k=i;k<m;k++) a[k,i] *= scale;
						}
					}
					w[i]=scale *g;
					g=s=scale=0.0;
					if (i+1 <= m && i+1 != n) 
					{
						for (k=l-1;k<n;k++) scale += Math.Abs(a[i,k]);
						if (scale != 0.0) 
						{
							for (k=l-1;k<n;k++) 
							{
								a[i,k] /= scale;
								s += a[i,k]*a[i,k];
							}
							f=a[i,l-1];
							g = -Sign(Math.Sqrt(s),f);
							h=f*g-s;
							a[i,l-1]=f-g;
							for (k=l-1;k<n;k++) rv1[k]=a[i,k]/h;
							for (j=l-1;j<m;j++) 
							{
								for (s=0.0,k=l-1;k<n;k++) s += a[j,k]*a[i,k];
								for (k=l-1;k<n;k++) a[j,k] += s*rv1[k];
							}
							for (k=l-1;k<n;k++) a[i,k] *= scale;
						}
					}
					anorm=Math.Max(anorm,(Math.Abs(w[i])+Math.Abs(rv1[i])));
				}
				for (i=n-1;i>=0;i--) 
				{
					if (i < n-1) 
					{
						if (g != 0.0) 
						{
							for (j=l;j<n;j++)
								v[j,i]=(a[i,j]/a[i,l])/g;
							for (j=l;j<n;j++) 
							{
								for (s=0.0,k=l;k<n;k++) s += a[i,k]*v[k,j];
								for (k=l;k<n;k++) v[k,j] += s*v[k,i];
							}
						}
						for (j=l;j<n;j++) v[i,j]=v[j,i]=0.0;
					}
					v[i,i]=1.0;
					g=rv1[i];
					l=i;
				}
				for (i=Math.Min(m,n)-1;i>=0;i--) 
				{
					l=i+1;
					g=w[i];
					for (j=l;j<n;j++) a[i,j]=0.0;
					if (g != 0.0) 
					{
						g=1.0/g;
						for (j=l;j<n;j++) 
						{
							for (s=0.0,k=l;k<m;k++) s += a[k,i]*a[k,j];
							f=(s/a[i,i])*g;
							for (k=i;k<m;k++) a[k,j] += f*a[k,i];
						}
						for (j=i;j<m;j++) a[j,i] *= g;
					} 
					else for (j=i;j<m;j++) a[j,i]=0.0;
					++a[i,i];
				}
				for (k=n-1;k>=0;k--) 
				{
					for (its=0;its<30;its++) 
					{
						flag=true;
						for (l=k;l>=0;l--) 
						{
							nm=l-1;
							if (Math.Abs(rv1[l])+anorm == anorm) 
							{
								flag=false;
								break;
							}
							if (Math.Abs(w[nm])+anorm == anorm) break;
						}
						if (flag) 
						{
							c=0.0;
							s=1.0;
							for (i=l;i<k+1;i++) 
							{
								f=s*rv1[i];
								rv1[i]=c*rv1[i];
								if (Math.Abs(f)+anorm == anorm) break;
								g=w[i];
								h=Pythag.Run(f,g);
								w[i]=h;
								h=1.0/h;
								c=g*h;
								s = -f*h;
								for (j=0;j<m;j++) 
								{
									y=a[j,nm];
									z=a[j,i];
									a[j,nm]=y*c+z*s;
									a[j,i]=z*c-y*s;
								}
							}
						}
						z=w[k];
						if (l == k) 
						{
							if (z < 0.0) 
							{
								w[k] = -z;
								for (j=0;j<n;j++) v[j,k] = -v[j,k];
							}
							break;
						}
						if (its == 29)
							throw new Exception("No convergence in 30 iterations of Svd.Decomp");
						x=w[l];
						nm=k-1;
						y=w[nm];
						g=rv1[nm];
						h=rv1[k];
						f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
						g=Pythag.Run(f,1.0);
						f=((x-z)*(x+z)+h*((y/(f+Sign(g,f)))-h))/x;
						c=s=1.0;
						for (j=l;j<=nm;j++) 
						{
							i=j+1;
							g=rv1[i];
							y=w[i];
							h=s*g;
							g=c*g;
							z=Pythag.Run(f,h);
							rv1[j]=z;
							c=f/z;
							s=h/z;
							f=x*c+g*s;
							g=g*c-x*s;
							h=y*s;
							y *= c;
							for (jj=0;jj<n;jj++) 
							{
								x=v[jj,j];
								z=v[jj,i];
								v[jj,j]=x*c+z*s;
								v[jj,i]=z*c-x*s;
							}
							z=Pythag.Run(f,h);
							w[j]=z;
							if(z!=0.0) 
							{
								z=1.0/z;
								c=f*z;
								s=h*z;
							}
							f=c*g+s*y;
							x=c*y-s*g;
							for (jj=0;jj<m;jj++) 
							{
								y=a[jj,j];
								z=a[jj,i];
								a[jj,j]=y*c+z*s;
								a[jj,i]=z*c-y*s;
							}
						}
						rv1[l]=0.0;
						rv1[k]=f;
						w[k]=x;
					}
				}
			}
		}
		public class SvdFloat
		{
			private static int MAX_ITERATIONS = 30;
			private static float anorm,c,f,g,h,s,scale,x,y,z;
			private static float[] rv1, tmp;
			private static int i,its,j,jj,k,l=0,nm=0;
			private static bool flag;
			/// <summary>
			/// Performs singular value back-substitution into x (A=U*diag(w)*V') after Svdcmp.
			/// </summary>
			/// <param name="u">Input matrix, the "a" argument after Svdcmp was run on it.</param>
			/// <param name="m">The number of rows in a.</param>
			/// <param name="n">The number of columns in a.</param>
			/// <param name="w">Output vector from Svdcmp.</param>
			/// <param name="v">Output vector from Svdcmp.</param>
			/// <param name="b">Right hand side input.</param>
			/// <param name="x">Solution vector output.</param>
			public static void Backsub(float[,] u, int m, int n, float[] w, float[,] v, float[] b, float[] x)
			{
				if(tmp==null || tmp.Length<n)
					tmp=new float[n];
				for (j=0;j<n;j++) 
				{
					s=0.0f;
					if (w[j] != 0.0) 
					{
						for (i=0;i<m;i++) s += u[i,j]*b[i];
						s /= w[j];
					}
					tmp[j]=s;
				}
				for (j=0;j<n;j++) 
				{
					s=0.0f;
					for (jj=0;jj<n;jj++) s += v[j,jj]*tmp[jj];
					x[j]=s;
				}
			}
			/// <summary>
			/// Performs the Singular Value Decomposition of "a" (matrix A), which is replaced by matrix U
			/// on output (for input into Svdbks).  w is the diagonal component of W.  v is V.  A=UWV'.
			/// </summary>
			/// <param name="a">The matrix to be decomposed on input.  On output, the matrix U.</param>
			/// <param name="m">int  The number of rows in a.</param>
			/// <param name="n">int  The number of columns in a.</param>
			/// <param name="w"></param>
			/// <param name="v"></param>
			public static void Decomp(float[,] a, int m, int n, float[] w, float[,] v)
			{	
				if(rv1==null || rv1.Length<n)
					rv1=new float[n];
				g=scale=anorm=0.0f;
				for (i=0;i<n;i++) 
				{
					l=i+2;
					rv1[i]=scale*g;
					g=s=scale=0.0f;
					if (i < m) 
					{
						for (k=i;k<m;k++) scale += Math.Abs(a[k,i]);
						if (scale != 0.0) 
						{
							for (k=i;k<m;k++) 
							{
								a[k,i] /= scale;
								s += a[k,i]*a[k,i];
							}
							f=a[i,i];
							g = -Sign((float)Math.Sqrt(s),f);
							h=f*g-s;
							a[i,i]=f-g;
							for (j=l-1;j<n;j++) 
							{
								for (s=0.0f,k=i;k<m;k++) s += a[k,i]*a[k,j];
								f=s/h;
								for (k=i;k<m;k++) a[k,j] += f*a[k,i];
							}
							for (k=i;k<m;k++) a[k,i] *= scale;
						}
					}
					w[i]=scale *g;
					g=s=scale=0.0f;
					if (i+1 <= m && i+1 != n) 
					{
						for (k=l-1;k<n;k++) scale += Math.Abs(a[i,k]);
						if (scale != 0.0) 
						{
							for (k=l-1;k<n;k++) 
							{
								a[i,k] /= scale;
								s += a[i,k]*a[i,k];
							}
							f=a[i,l-1];
							g = -Sign((float)Math.Sqrt(s),f);
							h=f*g-s;
							a[i,l-1]=f-g;
							for (k=l-1;k<n;k++) rv1[k]=a[i,k]/h;
							for (j=l-1;j<m;j++) 
							{
								for (s=0.0f,k=l-1;k<n;k++) s += a[j,k]*a[i,k];
								for (k=l-1;k<n;k++) a[j,k] += s*rv1[k];
							}
							for (k=l-1;k<n;k++) a[i,k] *= scale;
						}
					}
					anorm=Math.Max(anorm,(Math.Abs(w[i])+Math.Abs(rv1[i])));
				}
				for (i=n-1;i>=0;i--) 
				{
					if (i < n-1) 
					{
						if (g != 0.0) 
						{
							for (j=l;j<n;j++)
								v[j,i]=(a[i,j]/a[i,l])/g;
							for (j=l;j<n;j++) 
							{
								for (s=0.0f,k=l;k<n;k++) s += a[i,k]*v[k,j];
								for (k=l;k<n;k++) v[k,j] += s*v[k,i];
							}
						}
						for (j=l;j<n;j++) v[i,j]=v[j,i]=0.0f;
					}
					v[i,i]=1.0f;
					g=rv1[i];
					l=i;
				}
				for (i=Math.Min(m,n)-1;i>=0;i--) 
				{
					l=i+1;
					g=w[i];
					for (j=l;j<n;j++) a[i,j]=0.0f;
					if (g != 0.0) 
					{
						g=1.0f/g;
						for (j=l;j<n;j++) 
						{
							for (s=0.0f,k=l;k<m;k++) s += a[k,i]*a[k,j];
							f=(s/a[i,i])*g;
							for (k=i;k<m;k++) a[k,j] += f*a[k,i];
						}
						for (j=i;j<m;j++) a[j,i] *= g;
					} 
					else for (j=i;j<m;j++) a[j,i]=0.0f;
					++a[i,i];
				}
				for (k=n-1;k>=0;k--) 
				{
					for (its=0;its<MAX_ITERATIONS;its++) 
					{
						flag=true;
						for (l=k;l>=0;l--) 
						{
							nm=l-1;
							if (Math.Abs(rv1[l])+anorm == anorm) 
							{
								flag=false;
								break;
							}
							if (Math.Abs(w[nm])+anorm == anorm) break;
						}
						if (flag) 
						{
							c=0.0f;
							s=1.0f;
							for (i=l;i<k+1;i++) 
							{
								f=s*rv1[i];
								rv1[i]=c*rv1[i];
								if (Math.Abs(f)+anorm == anorm) break;
								g=w[i];
								h=Pythag.Run(f,g);
								w[i]=h;
								h=1.0f/h;
								c=g*h;
								s = -f*h;
								for (j=0;j<m;j++) 
								{
									y=a[j,nm];
									z=a[j,i];
									a[j,nm]=y*c+z*s;
									a[j,i]=z*c-y*s;
								}
							}
						}
						z=w[k];
						if (l == k) 
						{
							if (z < 0.0) 
							{
								w[k] = -z;
								for (j=0;j<n;j++) v[j,k] = -v[j,k];
							}
							break;
						}
						if (its == MAX_ITERATIONS-1)
							throw new Exception("No convergence in "+Convert.ToString(MAX_ITERATIONS)+" iterations of Svd.Decomp");
						x=w[l];
						nm=k-1;
						y=w[nm];
						g=rv1[nm];
						h=rv1[k];
						f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0f*h*y);
						g=Pythag.Run(f,1.0f);
						f=((x-z)*(x+z)+h*((y/(f+Sign(g,f)))-h))/x;
						c=s=1.0f;
						for (j=l;j<=nm;j++) 
						{
							i=j+1;
							g=rv1[i];
							y=w[i];
							h=s*g;
							g=c*g;
							z=Pythag.Run(f,h);
							rv1[j]=z;
							c=f/z;
							s=h/z;
							f=x*c+g*s;
							g=g*c-x*s;
							h=y*s;
							y *= c;
							for (jj=0;jj<n;jj++) 
							{
								x=v[jj,j];
								z=v[jj,i];
								v[jj,j]=x*c+z*s;
								v[jj,i]=z*c-x*s;
							}
							z=Pythag.Run(f,h);
							w[j]=z;
							if(z!=0.0) 
							{
								z=1.0f/z;
								c=f*z;
								s=h*z;
							}
							f=c*g+s*y;
							x=c*y-s*g;
							for (jj=0;jj<m;jj++) 
							{
								y=a[jj,j];
								z=a[jj,i];
								a[jj,j]=y*c+z*s;
								a[jj,i]=z*c-y*s;
							}
						}
						rv1[l]=0.0f;
						rv1[k]=f;
						w[k]=x;
					}
				}
			}
		}
		public class Pythag
		{
			static double absa,absb;
			/// <summary>
			/// Given the sides of a triangle, computes the hypotenuse "c" without the possibility
			/// of generating an exception.
			/// </summary>
			/// <param name="a">One side of a triangle.</param>
			/// <param name="b">Another side of a triangle.</param>
			/// <returns></returns>
			public static double Run(double a, double b)
			{
				absa=Math.Abs(a);
				absb=Math.Abs(b);
				if (absa > absb) return absa*Math.Sqrt(1.0+Math.Pow(absb/absa, 2.0));
				else return (absb == 0.0 ? 0.0 : absb*Math.Sqrt(1.0+Math.Pow(absa/absb, 2.0)));
			}
			/// <summary>
			/// Given the sides of a triangle, computes the hypotenuse "c" without the possibility
			/// of generating an exception.
			/// </summary>
			/// <param name="a">One side of a triangle.</param>
			/// <param name="b">Another side of a triangle.</param>
			/// <returns></returns>
			public static float Run(float a, float b)
			{
				absa=(double)Math.Abs(a);
				absb=(double)Math.Abs(b);
				if (absa > absb) return (float)(absa*Math.Sqrt(1.0+Math.Pow(absb/absa, 2.0)));
				else return (absb == 0.0 ? 0.0f : (float)(absb*Math.Sqrt(1.0+Math.Pow(absa/absb, 2.0))));
			}
		}
		/// <summary>
		/// If b is positive, then return a.  Otherwise, return -a.
		/// </summary>
		/// <param name="a"></param>
		/// <param name="b"></param>
		/// <returns></returns>
		public static double Sign(double a, double b)
		{
			if(b>=0.0) return a;
			return -a;
		}
		/// <summary>
		/// If b is positive, then return a.  Otherwise, return -a.
		/// </summary>
		/// <param name="a"></param>
		/// <param name="b"></param>
		/// <returns></returns>
		public static float Sign(float a, float b)
		{
			if(b>=0.0f) return a;
			return -a;
		}
		/// <summary>
		/// Swap functionality;
		/// </summary>
		public class Swap
		{
			public static double dtemp;
			public static float ftemp;
			public static int itemp;
			public static void Run(ref double a, ref double b)
			{
				dtemp = a;
				a = b;
				b = dtemp;
			}
			public static void Run(ref int a, ref int b)
			{
				itemp = a;
				a = b;
				b = itemp;
			}
			public static void Run(ref float a, ref float b)
			{
				ftemp = a;
				a = b;
				b = ftemp;
			}
		}
	}
}